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Plexcitons are polaritonic modes that result from the strong coupling between excitons and plasmons. We 
consider plexcitons emerging from the interaction of excitons in an organic molecular layer with surface 
plasmons in a metallic film. We predict the emergence of Dirac cones in the two-dimensional bandstructure 
of plexcitons due to the inherent alignment of the excitonic transitions in the organic layer. These Dirac 
cones may open up in energy by simultaneously interfacing the metal with a magneto-optical layer and 
subjecting the whole system to a perpendicular magnetic field. The resulting energy gap becomes populated 
with topologically protected one-way modes which travel at the interface of this plexcitonic system. Our 
theoretical proposal suggests that plexcitons are a convenient and simple platform for the exploration of 
exotic phases of matter as well as of novel ways to direct energy flow at the nanoscale. 


When UV-visible light is absorbed by an organic molec¬ 
ular aggregate, it promotes molecules from their ground 
to their excited electronic states. The resulting excitations, 
known as excitons, can migrate between molecules via a 
mixture of coherent and incoherent processes [ll]. Under¬ 
standing and controlling how this migration of energy oc¬ 
curs is a fundamental problem of chemistry and physics of 
condensed phases. Furthermore, it is also a technological 
problem which is relevant to the development of efficient 
organic solar cells and light-emitting devices as well as all- 
optical circuitry @1 . Many strategies to enhance the motion 
of excitons exist, a particularly interesting one being where 
they couple to surface plasmons (SPs) i^- In such strat¬ 
egy, the spatial coherence of plasmons assists the transport 
of an exciton across lengthscales that are dozens of times 
larger than regular exciton diffusion lengths. When the cou¬ 
pling is strong, meaning that the energy exchange between 
the exciton and plasmon is faster than the respective decay 
times la.laIg-Ulll . plexcitons (a class of polaritons) emerge 
d,[l^, and energy can migrate ballistically over the coher¬ 
ence length of the plasmon. Besides their usefulness in en¬ 
ergy transport, organic plexcitons are promised to be an ex¬ 
citing room-temperature "laboratory” for the study of light- 
matter and many-body interactions at the nanoscale |Tl|] . 
In this letter, we propose novel plexcitonic phenomenona 
which should be readUy realizable with current experimen¬ 
tal capabilities: Dirac cones and topologically nontrivial 
plexcitons which travel along preferred directions at the 
edge of an organic layer. 

Topologically nontrivial states of matter have been a topic 
of great interest in condensed matter physics owing to the 
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discoveries of the Quantum HaU Effect [HI]. and more re¬ 
cently, of topological insulators (I3.[3]- The systems sup¬ 
porting these states are characterized by topological invari¬ 
ants [ij], integer numbers that remain unchanged by weak 
perturbations. Physically, a nontrivial topological invariant 
signals the presence of one-way edge modes that are im¬ 
mune against moderate amount of disorder. Even though 
these phenomena were first conceptualized for fermions in 
solids, they have been successfully generalized to bosonic 
systems including photons in waveguides EM] , ring res- 
onator arrays (13. ultracold atoms in optical lattices [13. 
and classical electric circuits [^,[H|]- Furthermore, we have 
recently proposed an excitonic system consisting of a two- 
dimensional porphyrin film which becomes topologically 
nontrivial in the presence of a magnetic field . A chal¬ 
lenging feature of that proposal is the requirement of large 
magnetic fields (>10 T) and cryogenic temperatures to pre¬ 
serve exciton coherence. Even though we do not discour¬ 
age the experimental implementation of the latter, we con¬ 
sider a conceptually different platform which, by using plex¬ 
citons, avoids the use of large magnetic fields and, under 
appropriate circumstances, may work at room temperature. 
In the last year, Dirac and topological polaritons have been 
proposed in other contexts, such as optomechanical ar¬ 
rays and inorganic materials in optical cavities. All of these 
works share a common goal to ours, which is the design 
of exotic modes in strongly coupled light-matter systems. 
However, there are substantial qualitative and quantitative 
differences arising from the choices of material (organic 
exciton vs inorganic exciton Esi - lisll or mechanical mode 
[ 13113 ]) and electromagnetic (SP vs microcavity (I^ - Esll or 
photonic crystal [23,113]) excitations. Hence, the physics in¬ 
volved in our plexciton system contrasts with the other pro¬ 
posals in terms of the energy and lengthscales involved in 
the excitations, the magnitude of the couplings, the genera- 
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Figure 1. Plexciton setup. It consists of a plasmonic metal film, a 
dielectric spacer, and an organic layer. The latter is taken to be a 
monoclinic superlattice which makes an angle 6 with respect to 
the X axis, and is further characterized by angles P,J,S as well as 
distances between nanopillars A;, and Ay. Each element in the 
superlattice is a nanopiUar of organic emitters of dimensions W; 
along each axis. When the density of emitters is big enough, the 
coupling between the excitons in the organic layer and the surface 
plasmons (SPs) in the metal becomes larger than their linewidths, 
giving rise to polaritonic eigenmodes that are superpositions of ex¬ 
citons and plasmons, or more succintly, plexcitons. In this arti¬ 
cle, we shall also consider the case where the dielectric spacer is a 
magneto-optical (MO) material. 


tion of nontrivial topology, and the experimental conditions 
for its realization. Organic excitons differ from their inor¬ 
ganic counterparts in that they have large binding energies 
and are associated with large transition dipole moments. 
SP electromagnetic fields are strongly confined compared 
to their microcavity counterparts. The combination of all 
these properties in the organic plexciton context gives rise 
to strong light-matter interactions even at room tempera¬ 
ture and in an “open cavity” setup [H. 

The setup of interest is depicted in Fig. [H It consists of 
three layers, from bottom to top: a plasmonic metal mod- 

0)^ 

elled with a Drude permittivity = foo “ Coo ~ 4, 

cjp ~ 9 eV, representative parameters for Ag), an a = lOnm 
thick dielectric spacer (c^ ~ 1), and an organic layer (Corg ~ 
1). The spacer is placed to avoid quenching of the exci¬ 
tons by the surface plasmons (SPs) upon close contact, but 
as we shall see, we will also consider the case where it is a 
magneto-optical (MO) material. 

A quantum mechanical description of this setup is given 
by a Hamiltonian, 


H-Hexc + Hsp + Hexc-SP: ( 1 ) 

where each of the terms denotes the energetic contributions 
from the excitons in the organic layer, the SPs, and the cou¬ 
pling between them. More specifically (fi = 1), 


Hexc = Y.^nCT],CTn+ E + h.C.), (2a) 

n^n' 

= I^w(fc)afc«fc, (2b) 

k 

Hexc-SP^Yjc/knakt^W’^" + (2c) 

k,n 

Here, in order to obtain a specific shape of the exciton en¬ 
ergy dispersion (see below), we have taken the organic layer 
to be an oblique superlattice of organic nanopillars. The su¬ 
perlattice is defined by the distances A;, (horizontal) and 
Ay (vertical), as weU as angles p, y, 6, and 0 (Fig. [T^); 
the nanoprUars are in turn rectangular parallelepipeds of 
densely packed organic chromophores (assuming a van der 
Waals distance between chromophores of 0.3 nm, pnp =37 
chromophores/nm^) with volume V„p = FPyVPyWz (Fig.[T}D), 
obtained from growing a J-aggregate film . Given this 

preamble, cr],((T„) and a^(flfc) label the creation (annihila¬ 
tion) operators for the collective exciton at the n -th nanopil¬ 
lar and the fc-th SP mode, respectively, where n and k are 
(two-dimensional) in-plane vectors denoting a position and 
a wavevector, respectively. J-aggregation of chromophores 
results in a collective transition dipole at an excitation 
energy w„, while the dispersion energy of the fc-th SP mode 
is denoted wik). Dipolar interactions J„„i couple the vari¬ 
ous nanopillars. The coupling between the exciton and the 
SP depends on the average in-plane location r„ of the n-th 
nanopillar, and is also dipolar in nature. 


Jkn = J --BCfc)- (3) 

y 2eo SLk 

Here, Lfc denotes a vertical (z-direction) mode-length of the 
SP which guarantees that the total energy of a SP prepared at 
the fcth mode is quantized at the energy w(fc), E{k) is an ap¬ 
propriately scaled electric field of the corresponding mode, 
and g-“orgo(fc)z(fc) yields a mean-field average of the interac¬ 
tion of the evanescent SP field (with decay constant aorgo in 
the organic layer) over the chromophores at different verti¬ 
cal positions of the nanopillar; it optimizes the interaction 
such that one may assume the nanopillar is a point-dipole 
located at the mean height z{k). The latter average renders 
the originally 3D system into an effectively 2D one. Detailed 
derivations of Eqs. l|2) and 0 are available in the Supple¬ 
mentary Information (SI)-III. 

Assuming perfect periodicity of the superlattice (w„ = 
w, - fi) and only nearest and next-nearest neighbor 
(NN and NNN) dipolar interactions, we can re-express Hgxc 
(Eq. (2al ) in terms of k modes. As explained in SI-IIIA, 
it is possible to approximate Hexc up to 0(|fcp) as arising 
from an ejfective simplified rectangular (rather than mon¬ 
oclinic) lattice aligned along the x, y axes, and with NN 
interactions only. This is a reasonable thing to do, as the 
topological effects we are interested arise at relatively long 
wavelengths. Within this approximation, we may construct 
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Fourier modes at - , ^ Ft./ al.e (here, Ni is the 

k' 

effective number of nanopillars along the z-th direction) 
and rewriting Eq. (1) in reciprocal space, we obtain H - 
LkHk, 

JJfe = ^<Jeji:c,fcO'^crfc + w(fc)aj[.flfc+ [^(fc)flfcO’^ + h.c.], (4) 

= Hexc,k =Hsp,k =Hexc-SP,k 

where, 


<^exc,k = (^eff + 2JxCOS{kxAx) + 2/y COSffcyAy), (5a) 


where /; and A; are effective NN hoppings and spac- 
ings along the z-th axis, respectively, Wg// is the effective 
nanopillar site energy, and we have taken ^ for all n. 
Eq. I5bt denotes the fc-dependent coupling between a plas- 
monic mode and a collective exciton state throughout the 
organic layer. It features the interaction between nanopU- 
lars and the plasmonic mode, which is proportional to the 
square root of the total number of nanopillars yJNxNy co¬ 
herently coupled to the plasmon (compare with Eq. (^) Isol - 

m. 

For our simulation, we choose the length parameters 
Ah = 100, Ay = 88, Wx = 7.5, Wy = 50, = 70 nm, and 

take p - 13.1° (see SI-llLA for explanation of choice of pa¬ 
rameters). Denoting the transition dipole fi - fifi, we es¬ 
timate p = \/Nnp X lOD = 9855D, where Nnp = PnpVnp is 
the number of chromophores in the nanopillar; p, is the in¬ 
plane unit vector making an angle of a = 223°) with respect 
to X, that is, p - cosax -i- sinay. Choosing the simulation 
values for H^xc in Eq- tSal to be Aj,; = Ay = 50 nm, we get the 
effective parameters - 2.15eV, Jx = 362meV, and Jy - 
107meV. Fig. |2^ shows superimposed dispersion curves 
for Hexc,k and Hspk independently (taking Hexc-SP = 0). 
The superlattice has been intentionally constructed to ob¬ 
tain Jx, Jy > 0 yielding a “dome”-like dispersion for Hgxc.k, 
i.e. it features a maximum at fc = 0, behaving as a 2D “IT- 
aggregate” a . As for Hspk, its dispersion has the shape of 
a rotationally symmetric “fountain,” and is nothing more 
than the 2D rendering of the standard ID textbook result 
fl^l . featuring a linear dependence of the energy at short 
wavevectors and a plateau at large ones, indicating excita¬ 
tions that are qualitatively closer to light or to charge oscil¬ 
lations in the metal, respectively. 

Fig. 1^ shows the two-plexciton-branch bandstructure 
arising from the diagonalization of Eq. JT) . We notice that 
anticrossing gaps are noticeably opened in the vicinity of 
where the dispersion curves for Hexc,k and H^pk used to 
cross in Fig. |2^. This is a signature of exciton-SP coupling 
Hgxc.SP.k- Given a fixed wavevector direction, whenever 
these anticrossings occur, the lower-plexciton (LP) branch 
starts off as being mostly SP at short k values, but paramet¬ 
rically morphs into mosdy exciton at large ones; the oppo¬ 
site happens with the upper-plexciton (UP) branch. ITow- 
ever, probably the most striking feature of Fig. |2 )d is the 


appearance of two Dirac cones (see dashed circles) at crit¬ 
ical wavevectors k* in which anticrossings do not happen. 
Their onset coincides with the directions at which k is or¬ 
thogonal to fi. Their physical origin is explained in Fig. |2^ 
which, in its top panel, shows the in-plane electric field for 
the k-th SP mode, E±[k) = E[k) - E{k) ■ zz and is purely 
parallel to k, E±{k) - Ekik) (blue vector field). If all the 
dipoles in the organic layer are aligned (in-plane) along /x, 
their projection onto the SP electric field, which gives rise to 
the exciton-SP coupling (see Eq. (5b)), will wax and wane as 
a function of the azimuthal angle q) between the fixed dipole 
and the varying SP wavevector according to fi-E{k)(x coscp. 
Clearly, this projection will vanish if k happens to be orthog¬ 
onal to fi, that is, at the special angles (p - f,^, so that 
any degeneracy between the exciton and the SP modes will 
remain unlifted along these directions. From this physical 
picture, we can extract the two essential ingredients for the 
emergence of the plexciton Dirac cones. First, the dipoles 
need to be aligned to create an anisotropic exciton-SP cou¬ 
pling as a function of (p. Second, this alignment needs to be 
horizontal, as a vertical component of the dipole will cou¬ 
ple to the vertical component of the electric field £(fc) ■ zz, 
and this coupling, unlike its horizontal counterpart, does 
not vanish for any qi. It is important to note that neither 
of these requirements requires the use of the superlattice, 
which will be exploited for a different purpose (the engi¬ 
neering of topological edge modes, as explained in the next 
paragraph). Therefore, a standard organic molecular crys¬ 
tal with aligned transition dipoles lying on the horizontal 
xy plane will suffice. These plexciton Dirac cones which, to 
our knowledge, have not been reported in the past, should 
be easily detectable by collecting the reflected light spectra 
upon excitation of the plexcitonic system in a grating, Otto, 
or Kretschmann configurations [I^ , by systematically scan¬ 
ning across |fc| and qt values. For a general {\k\,q), the spec¬ 
trum should consist of two dips as a function of dispersed 
energy, each associated with the corresponding eigenener- 
gies of the LP and UP. Importantly, however, the two dips 
merge at the Dirac cones. In a standard plexciton disper¬ 
sion measurement, one only scans across |A:|. Since we are 
interested in a two-dimensional dispersion, the scan must 
also be performed across q. 

Having elucidated the mechanism for the formation of 
plexciton Dirac cones, we proceed to entertain a more am¬ 
bitious goal. We aim to engineer topologically protected 
plexcitons by opening the Dirac cones using a time-reversal 
symmetry breaking (TRSB) perturbation [ijl. To accom¬ 
plish this, we now assume that the dielectric spacer has 
magneto-optical (MO) properties; that is, upon application 
of a perpendicular magnetic field, its permittivity becomes 
anisotropic, mo, 


£ MO- 


£d ig 0 
-ig £d 0 
0 0 ea 


( 6 ) 


Materials associated with this dielectric tensor exhibit Fara¬ 
day effect and are of great interest in the fabrication of 
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Figure 2. Bulk plexciton properties. Dispersion relations: (a) for SP and exciton (organic layer) modes independently, (b) when they couple 
in the absence of the MO effect, yielding lower (LP) and upper (UP) plexciton branches which feature two Dirac cones (dashed hlack 
circles), and (c) when they couple in the presence of the MO effect (g = 0.3), lifting the Dirac cones, (d) Berry curvature associated with 
the LP in (c). (e) Physical mechanism for the appearance of plexciton Dirac points: (top) collective (in-plane) exciton transition dipole 
for a nanopillar; (bottom) magnitude of the electric field of magneto-SP modes as a function of wavevector. In the absence of the MO 
effect, only the wavevector-parallel components ijj.(blue) are present. Thus, the nanopillars experience no coupling with modes whose 
wavevectors are perpendicular to the transition dipole. Along these directions, degeneracies between the SP and the exciton modes are 
not lifted, yielding two plexciton Dirac points. Nonzero tangential components Eq (red) emerge upon inclusion of the MO effect, lifting 
these degeneracies. 


optical isolators . For the near-IR and visible, yttrium 
iron garnets (YaFesOiz, YIG) substituted with Bi (BiYlG), Ge 
(CeYIG) or other rare earths provide high Faraday rotation 
with low optical absorption [ula ; alternatively, one may use 
an MO active Co-alloy film or multilayer with an in¬ 
sulating layer to avoid quenching of the excitons. In this 
article, we are interested in the new SP modes, denoted 
as magneto-SP modes, arising at the interface of the plas- 
monic metal and the MO layer. The solution for the latter 
is highly non-trivial, and we refer the reader to our pertur¬ 
bative solution in the SI-1,11, which builds a perturbation 
theory based upon an initial calculation due to Chiu and 
Quinn Q. Fig. 1^ shows essentially the same calculation 
as Fig. except for the inclusion of the MO effect. In¬ 
terestingly, we notice that the Dirac cones have been lifted. 
A physical understanding of the latter phenomenon can be 
obtained by appealing to Fig. ^ again. Within our pertur¬ 
bation theory, the magneto-SP modes differ from their orig¬ 
inal SP counterparts in that there are additional tangential 
components (red) to the electric field. The clockwise vortex 
vector field is a signature of TRSB; it becomes counterclock¬ 
wise upon change of direction of the magnetic field. This 
tangential electric field is the sole responsible for opening 
the Dirac cones at the critical angles q) - f, where the 
original field (blue) ceased to couple to the excitons. Hence, 
we have concocted a situation where anticrossings occur for 
all azimuthal angles cp. To characterize the topology of the 
resulting bandstructure, we numerically compute the Berry 


curvature for each plexciton branch ; we show that of 
the LP in Fig. |2ll. Its integral with respect to the BriUouin 
zone is the so-caUed Ghern number C, an integer which, if 
nonzero, signals a topologically nontrivial phase. Fig. 
clearly shows that this integral is non-vanishing, and in fact, 
adds up to C = -1 (by the sum rule of Ghern numbers, the 
upper branch necessarily has C = 1). Intuitively, it is also 
clear that most of the nontrivial topology, and hence, Berry 
curvature, is concentrated in the vicinity of what used to 
be the Dirac cones. In passing, we note that considerable 
attention has recently been given to magneto-SPs where 
the magnetic field is applied parallel (instead of perpendic¬ 
ular) to the metal film itself, yielding dispersion relations 
which are nonreciprocal 0.’ Guriously, this arrangement 
does not give us the topological states we are looking for, 
although it might be intriguing to explore the connection 
between these magneto-SPs and the ones exploited in our 
present work, arising from a perpendicular magnetic field. 

So far, all the described calculations have been carried 
out in the bulk. By virtue of the bulk-boundary correspon¬ 
dence Q , we expect topologically protected one-way edge 
modes associated with this setup. In order to compute 
them, it is convenient to keep periodic boundary conditions 
for the magneto-SP modes, yet consider two domains of ex¬ 
citons on top (Fig. [ 3 ^), one with (in-plane) dipoles pointing 
along fi (red dipoles) and the other one with vertical dipoles 
along z (blue dipoles). Pictorially, this setup resembles a 
"donut with two icings,’’ where the donut is the metal with 
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toroidal geometry, and the two icings are the domains of ex- 
citons separated hy two interfaces located at y = + -^ and 
y = 0, where L; is the total width of the simulated sample 
along i (in our calculations, we take Lx - Ly - & fim). The 
Chern numbers associated with the hulk LP branch of each 
domain are C = -1 and C = 0, respectively. Hence, the plex- 
citons for the blue domain are topologically trivial. This can 
be understood by recalling that no plexciton Dirac points 
occur when dipoles are vertically aligned, regardless of the 
MO effect. In the limit of no disorder along the x-direction, 
kx is still a good quantum number, and Fig. shows the 
corresponding plexciton dispersion relation. This band- 
structure is essentially a projection of the gapped 2D bulk 
bandstructures of both domains of plexcitons onto one axis 
kx with additional states spanning the topological gap be¬ 
tween the LP and UP branches. Inspection of the nature of 
these mid-gap states reveals that they have substantial exci- 
ton and magneto-SP character, and that they are precisely 
the edge states we are searching for: one band has posi¬ 
tive (negative) dispersion and is localized along y = 0 (y = 
+ -^). Thus, by preparing a plexciton wavepacket localized 
along one of the interfaces, and making sure it is composed 
of energy states within this topological gap, one ensures 
that transport occurs robustly without much probability of 
backscattering. The reason being that backscattering re¬ 
quires coupling between counterpropagating modes which 
are separated by a distance -^, which is large compared 
to the width of the corresponding wavefunctions along y- 
Fig-llta,b) shows snapshots of the dynamics associated with 
these edge states. Panels (a,b) and (c,d) show a plexciton 
that starts localized at x = and x- - respectively, and 
tracks the one-way (to the left or to the right) nature of their 
motion within the femtosecond timescale. 

A useful ingredient guaranteeing the robustness of one¬ 
way transport of these edge modes is that they appear 
within a global gap, the latter of which is a consequence of 
our 2D H-aggregate superlattice design. It is easy to check 
that if /jc < 0 or Jy < 0, this global gap is not guaranteed 
anymore, and edge modes may become degenerate with 
bulk modes. Hence, these two types of modes could read¬ 
ily hybridize, yielding channels connecting one edge to the 
other, allowing for backscattering. Importantly however, the 
bandstructure remains topologically nontrivial even in the 
absence of such global gap, so even if perfect one-way trans¬ 
port is not observed in these cases, signatures of the latter 
may remain. Furthermore, we have shown that a superlat¬ 
tice with in-plane dipoles and NN and NNN couplings be¬ 
tween nanopillars gives rise to an effective 2D H-aggregate. 
It is unclear whether this superlattice is necessary for our 
goal, or whether a simple organic molecular crystal can 
yield a similar behavior once we take into account aU the 
dipolar contributions, from short to long range. These is¬ 
sues will be explored in future work. In the mean time, it suf¬ 
fices to note that, as a proof of concept, a global gap which 
hosts topologically protected edge states can be obtained by 
using an organic superlattice. 

Notice that the dispersion of the edge plexcitons is such 


that a subset falls within the light cone (w = ckx! \/e^ so 
far field excitation and detection of this fraction is possi¬ 
ble via interaction with the organic layer itself. The rest of 
the plexcitons can be probed using the already mentioned 
SP measurement techniques, by launching plexcitons excit¬ 
ing the metallic layer itself. Furthermore, the ballistic and 
one-way nature of these modes can in principle be eluci¬ 
dated using fluorescence microscopy |^. It is important 
to note that, owing to the topological nature of these states, 
perfect lattices are not required, so this phenomenon holds 
as long as orientational and site energy disorder induce per¬ 
turbations which are smaller than the topological anticross¬ 
ings. We tested these ideas by simulating lattices with disor¬ 
der in the site energies ((hgff ^eff + ^cbeff) as well as in 
the orientations of the dipoles (/i ^ cosA,^[cos(a; -i- Aa)x -i- 
sin(a -I- Aa)y] -i- sinA^z), where Aj are chosen to be Gaus¬ 
sian random variables centered at 0 and having disorder 
widths aj for each j = We//> <P- By systematically varying 

these widths independently and keeping track of the pres¬ 
ence of the one-way edge states, we noticed that the latter 
survive under large amounts of disorder, whose thresholds 
are approximately located at A(S^^^ ~ 0.25 eV, A^ ~ 30”, and 

To clearly illustrate the nature of the topologically pro¬ 
tected edge states, we have taken c,/ = 1, g = 0.3, yield¬ 
ing a minimum gap between plexciton branches (at the 
wavevectors k* of the original Dirac points) of 0.48 eV The 
crossing of the SP and exciton dispersion curves happens 
at 3.1 eV Given typical linewidths associated with the vari¬ 
ous dissipative mechanisms at room temperature {jexc.rel ~ 
5meV, Jexc,deph ~ 40meV, jsp.rel ~ lOmeV, where rel and 
deph stand for relaxation and dephasing), we anticipate 
that exciton-magneto-SP couplings need to be at least 
|hfejcc-sp(fc*)| lOmeV in order for the topological plex¬ 
citon edge states to be meaningful {i.e., be in the regime 
of strong coupling [l^])- Since our perturbative theory is 
linear in g, this means that keeping all other parameters 
fixed, we require g > 0.03. Even though this value of MO ef¬ 
fect is very reasonable for BiYIG films, their polarizability is 
much larger than what we have used in our calculations; in 
fact, Erf ~ 6.25 and g ~ 0.1 when the perpendicular magnetic 

field is about 0.01 T [ll, yielding |Hejcc-Sp(fc*)| ~ lOmeV, 
which would render the topologic^ edge states difficult to 
detect, i.e., we need to reach the plexciton regime . A pos¬ 
sible solution to this problem is to consider novel MO gar¬ 
net compositions which maximize the g/e^ ratio. In SI- 
III-4, we summarize our understanding of the optimization 

of |Hej:c-Sp(fc*)| for the parameter space comprised by e^, 
g, and a. In our simulations, we have assumed that the 
organic-layer has been embedded into the MO layer, since 
the calculation with the proper three-layer setup yielded 
couplings |Hej:c-Sp(fc*) I that did not surpass the dissipative 
linewidths. However, as far as we are aware, the weakness 
of the MO effect obtained with MO garnets is not funda¬ 
mentally limited. Once we identify the physics that con¬ 
trols the strength of Hexc-spik *)\, we can optimize the MO 
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Figure 3. Topologically protected edge modes, (a) Simulation of 
edge modes where magneto-SPs are computed in the torus geom¬ 
etry. Two domains of organic layers are placed on top of it (just like 
"two icings” on a donut). In-plane (red) and out-of-plane (blue) 
transition dipoles yield topologically nontrivial and trivial plexci- 
tons, respectively. Topologically protected one-way plexcitons ap¬ 
pear at the interfaces (thick black arrows). Each interface features a 
different plexciton direction of motion, (b) ID dispersion relation 
(i)(kx) for the setup in (a). Bulk LPs (blue) and UPs (red) separated 
by edge modes (green) featuring positive and negative dispersions, 
respectively, and localized along each interface. 


layer by appealing to different engineering strategies to in¬ 
crease g and decrease (e.g. by using of MO garnet sphere 
arrays li^l, plasmonic/magnetic metal nanostructures I21I] , 
Ce substituted YIGs [3, Eu nanocrystals 13, etc.). 

In terms of the fabrication of the plexciton setup, we warn 
that the creation of BiYlG layers typically need high temper¬ 
ature and oxygen, which is incompatible with deposition on 
Ag or organic materials. Hence, the MO layer should first be 
deposited on garnet substrates such as GGG (Gd 3 Ga 50 i 2 ) 
(111), and subsequently floated off by dissolving or polish¬ 
ing the substrate. One should then transfer the film on a 
Ag-coated substrate and the organic layer may be deposited 
and patterned on garnet. 

To summarize, we have described the design of exotic 
plexcitons via a judicious choice of material and electro¬ 
magnetic excitation modes. We showed that Dirac cones 
and topologically protected edge states emerge from rel¬ 
atively simple hybrid organic/inorganic nanostructures. 
Even though we have not precisely identified an explicit 
MO material which fully satisfies our requirements, we be¬ 
lieve its design is within reach, and is the subject of our 
present investigations. It is also worth noting that the phys¬ 
ical origin of the described edges states is different from 
that of edge plasmons in disk geometries 13141], although 
the connections are intriguing. The possibility of directed 
migration of excitation energy at the nano- and mesoscale 
offers exciting prospects in light-harvesting and all-optical 
circuit architectures. Furthermore, given the recent ex¬ 


perimental discovery of nonlinear many-body effects such 
as Bose-Einstein condensation of organic cavity-polaritons 
l3 - [49ll and plexcitons jsl] at room temperature, the in¬ 
troduction of the novel features described in this letter en¬ 
riches the scope of these materials as a testbed for novel 
many-body quantum phenomena. 
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Supplementary Information for “Plexcitons: Dirac points and 

topological modes” 


In this Supplementary Information, we derive some results that are used in the main text of the article. We compute the 
electromagnetic fields associated with the magneto-surface-plasmon (magneto-SP] modes arising in two- and three-layer 
setups (Sec. [Jand Sec. [11) . Finally, we derive the effective 2D Hamiltonian for the plexciton system hy considering the dipolar 
couplings between chromophores in the organic layer and between the latter and the magneto-SPs (Sec. [111) . 


CONTENTS 


References 


s 


I. Magneto-SPs at a metal-dielectric interface 

A. Maxwell’s equations 

B. Permittivities 

C. Electromagnetic modes for each layer 

1. MO layer (z>0} 

2. Metal layer (z < 0) 

D. Matching the modes at the boundary (z = 0) 

E. Perturbation expansion on g 

1. Solving for Q 

2. Solving for 

3. Collecting the expressions for the fields 

F. Quantization and normalization of modes 

G. Perturbation expansion of Lfe 


S 

S 

_l 

1C 

1C 

U 

12 

13 

13 

15 

le 

17 

IS 


II. Magneto-SPs in a three-layer setup 

A. Electromagnetic modes for each layer 

1. Organic layer (z > a) 

2. MO layer (a > z > 0) 

3. Metal layer (z < 0) 

B. Matching the modes at the boundaries (z = 0 and z- a) 

C. Perturbation expansion on g 

1. Solving for t* and 

2. Collecting the expressions for the fields 

D. Quantization and normalization of modes 


ii 

2C 

21 

22 

22 

22 

24 

25 


III. Exciton-exciton and exciton-SP couplings 

A. Dipolar couplings between nanopillars 

B. Exciton-SP couplings 

1. Mean-field approximation (MEA) 

2. Beyond the MFA 

3. Comparison 

4. Representative coupling values 


26 

26 

2E 

32 


References 









z 



y 




■»x 


Figure SI. Metal-MO dielectric interface atz-0. We solve for the surface plasmon (SP) modes arising at the metal-dielectric interface, in 
particular, when the permittivity of the dielectric mo is anisotropic due to the application of a perpendicular external magnetic field. 
The SP modes in this case are referred to as magneto-SP modes. As a first approximation to the exciton-SP coupling, we assume that the 
organic layer is embedded within the dielectric spacer medium. 


I. MAGNETO-SPS AT A METAL-DIELECTRIC INTERFACE 
A. Maxwell’s equations 

Maxwell’s equations in arbitrary media read as follows. 


V-D = 0, (SI) 

V-.B = 0, (S2) 

VxE^-dtB, (S3) 

S/xH^dtD. (S4) 


We are interested in the case where the electric displacement D and the magnetic field H are related to the corresponding 
electric fields E and the magnetic inductions B via the linear constitutive relations D = cq'TE and H = where 

eo = 8.854 x and po = 1.257 x are the permittivity and permeability of vacuum, and and 'Ji are the 

corresponding scaling (unitless) tensors for the medium in question. By taking the curl of Eq. IS3t . we obtain the wave 
equation. 


V(V-.g)-V2E=-4-d2£, (S5) 

c2 

where we have introduced the free space speed of light c - . A completely analogous equation holds for H by taking 

the curl of Eq. tS4t , but the latter suffices for our purposes. 


B. Permittivities 


We are interested in (magneto-)SPs arising at the interface between a plasmonic metal (silver or gold) and a dielectric 
medium which is endowed with magneto-optical (MO) properties (see Fig. [^. The latter can be a mix of a magnetic oxide 
dissolved in a polymer. We assume that the magnetization of both layers is null, p = 1. 

The permittivity of the metal is taken to be isotropic ^ I , where 



0)2 -I- iwy 


£mito) ~ ^oo 


(S6) 
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is of the Drude form (the parameters for Ag (Au) are £^0 = 3.7(6.9), wp - 9.2(8.9) eV, and 7 = 0.01(0.07) eV Q]). Throughout 
this work, we will set 7 = 0 in order to keep the formalism simple. Physically, as long as the relevant energy scales of interest 
are larger than 7 (see main text), this is a good approximation. Since we are interest in the plexciton (strong exciton-SP 
coupling) regime IMJllll. this should he a good approximation to take. Otherwise, the quantization of the problem becomes 
much more complicated. 

For the MO layer, the permittivity is anisotropic: upon interaction with an external magnetic field in the perpendicular 
z-direction, it acquires the form. 


€ MO- 


Ed ig 0 
-ig Ed 0 

0 0 £d 


(S7) 


where the tensor has been written in Cartesian coordinates {x,y,z), and we take = 1 and g = 0.1. Here, the off-diagonal 
term is proportional to the Faraday rotation that a linearly polarized plane-wave electric field experiences as it passes 
through the material; g changes sign t^on change of magnetic field direction. Typically, MO magnetic oxides like Bismuth- 
and Yttrium-Iron Garnets (BIG, YIG) [JQl have permittivities of ~ 6 EUl , which imply a severe index mismatch with the 
metal and organic layers of interest. Hence, we are implicitly assuming that we have a Maxwell garnet blend with a low index 
polymer or aerogel [l^ at our disposition, which yields an effective e^-l. On the other hand, g ~ 0.1 is a reasonable param¬ 
eter for MO garnets under a magnetic field of 0.1 Tesla [ill. Just as with Cm, we ignore imaginary (absorptive) contributions 
to Erf. Subsec. IIIIB4ldis cusses other MO materials that could be used for the purposes of our study. 

Chiu and Quinn have solved a slightly different problem, namely, the magneto-SPs arising from a metal under a 
strong magnetic field coupled to a isotropic non-MO dielectric. In their study, the anisotropy arises in the metal permittivity 
rather than in the one corresponding to the dielectric. The resulting equations are, as expected, very similar, and one could 
translate their equations to our setup by some careful changes of variables. However, for clarity of presentation and in order 
to develop the three-layer calculation of Sec. (II), which is a generalization of the two-layer case, we shall outline the entire 
procedure here. Importantly, in doing so, we manage to go further than Chiu and Quinn and construct a perturbation theory 
in the small parameter g. This allows us to develop explicit expressions for the electromagnetic modes which, as far as we 
are aware, have not appeared in the literature before. 


C. Electromagnetic modes for each layer 

The problem is rotationaUy symmetric about the vertical z-direction, so it is convenient to adopt a cylindrical coordinate 
system. Let us search for SP modes labeled by k which propagate in-plane and decay along z (this is precisely the condition 
for SP modes). 


= (S 8 a) 

= (S 8 b) 

For a given direction of k, we shall write vectors in the right-handed cylindrical coordinate system spanned by the unit 
vectors zsuchthat A:x0(. = z. for instance, E - {Ek,Egi^,Ez), where J?,- = E-i (beware that we have defined the tangential 
direction of 0^ with respect to k and not to f). Physically, k and cj-wik) denote the in-plane (propagating) wavevector and 
frequency of the monochromatic wave, respectively, = r ■ fc is the projection of the position vector r = (r^, rg,z] along 
the k direction, and kz is the imaginary wavevector associated with the evanescent wave along the perpendicular direction. 
Inserting these modes into Eq. tSSi yields an anisotropic wave equation for E, 

jm~ jO^Ejki + kl — 0, (S9) 

Ijm 

where we have used the notation 0(A:) = corresponding to the free space wavevector. Here, the indices run 

through k,6,z and, formally, we may write k-{k]c,kg,kz)- {k, 0, fc^). Due to rotational symmetry, ”mo has the same form 
in the fc,0fc,z coordinates as Eq. (S7). 


1. MO layer (z > 0) 

Inserting the dielectric tensor associated with the MO layer (Eq. IS7l ) into Eq. (S9) yields a matrix equation Ml mo = 0 
which explicitly reads as. 
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Ed-El ^kl j^Q ig El ^kkz,MO 


1 


0 

-ig €d-El~^ik'^ + kl^^Q) 0 


Ee.MO 

- 

0 

kkz,MoEl~^ 0 £d-El~^k^ 


Ez.MO, 


0 


(SIO) 


At this point, we have chosen the arbitrary normalization condition Eic,mo - 1- The secular equation corresponding to Eq. 
(STO) is, 


where. 




+ C = 0, 


(Sll) 


IB = 

C = 




2 2 
£d 


(S12a) 

(S12h) 


The hi-quadratic Eq. (sn) yields solutions kz,MO - 


exponentially rising constants given hy. 


where are two different evanescent decay or 


“mo 



(S13) 


Note that these (in general, complex-valued) constants must have positive real part for e““Mo^ to decay or for e“Mo^ to 
rise, respectively. In the first two-layer setup we are considering, we will assume that the MO layer extends indefinitely for 
z > 0 so the field in this region must he a superposition of the two evanescent fields; exponentially rising fields wiU become 
important when we add an additional interface at z = a (see Sec. HD- The tangential and perpendicular components of the 
electric field (given E^^mo - 1) be obtained from Eq. ISIOI . 


F± = 

;■ Irrr- 

+ _ 

^z.MO- k2-Q2ed 


(S14a) 

(S14b) 


2. Metal layer (z < 0) 


In the metal, MmEm - 0 corresponds to, 


Em El 0 Q ^kkz,m 


1 


0 

0 £m-Er^{k^ + kl„i) 0 


^6,m 

- 

0 

kkz,mEl~^ 0 £m-El~^k^ 


J^z,m 


0 


where again we have chosen = 1- This leads to the secular equation which yields the exponentially decaying field for 
z < 0; by letting - -iocm, 


Eq. jSlSI reveals that 



(S16) 


ik 

Ez.m =- (S17) 

but does not inform us about the tangential component Eg^rn as the entry (Mlm)ee = £m “ ET^ik^ + m) - 0. This missing 
component will be deduced by matching the electromagnetic fields at the boundary z = 0. 
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D. Matching the modes at the boundary (z = 0) 

We are looking for magneto-SP modes labeled by a propagating wavevector k and frequency (w(fc), but, in general, different 
kz evanescent wavevectors for each layer, which we have denoted kr^^^ - ioL%[Q and kz,m - -iccm- The energy w{k) (and 
therefore O) is unknown. 

To summarize, for each pair [k,(i)), there are two possible modes in the MO layer associated with different decay constants 
oc%[Q (see Eq. IS13> , z > 0), 


~ ’ (S18a) 

^MO ~ 

-—[a-EgMO'~^^^z,MO~^MO’^^^0,MO^^^ ’ (S18b) 

where the magnetic induction in Eq. iS18bt has been deduced from Eq. IS18a) and Maxwell’s Eq. JSS). The vectors have 
been written in cylindrical coordinates and we have defined rj = jqgnce, the total fields in the MO layer read, 

EmO = ^MO^MO ^MO^MO’ 

Bmo - ^moBmo + ^mqBmo- 

where f- are coefficients to be determined. Similarly, for the metal layer (z < 0), 

Em — EmTjS 

= (l,E0,^,£z,m)^7e“"’^ (S20a) 

Bm — Bmf]S 

— {—(XmE^ m’~ ^^Ez m ^ ■ (S20b) 

( 1 ) 

Here, we keep the arbitrary normalization where Efc.m = 1. Later on, we shall fix this normalization via quantization of the 
energy of the modes (see|TF). Furthermore, it is also safe to arbitrarily assume E^^^ - 1 because Eqs. IS19al and IS19bt 
contain scaling coefficients f- which will be fixed by the boundary conditions at the metal-MO interface. 

We are ready to match the fields at the interface at z = 0. The in-plane electric field and the perpendicular electric displace¬ 
ment each need to be continuous across the boundary: Ei^mo - Ei,m for i - r,d, while eijEz,MO - £mEz,m- Furthermore, the 
magnetic field, and because 7? = 1, its induction, are aU continuous throughout, Bi MO - Bi^m - These constraints altogether 
read. 


(S19a) 

(S19b) 


kEz 


~ ^MO ^MO’ 

(S21a) 

~ ^MO^e.MO ^MqBq.MO’ 

(S21b) 

= ^dit^oEz,MO + ^MqBz.MqB 

(S21c) 

~ ~ ^MO^MoBb^mO ~ ^MO^MoBe,MO’ 

(S21d) 

- tMO^^Ez^MO - + ^MO^^Bz^mO ~ 

(S21e) 

- ^MqBb.MO ^MoBb.MO- 

(S21f) 


These constraints read similarly to the ones derived by Chiu and Quinn in [3] (soe their Eqs. (39)-(44)), except for the 
different coordinate conventions. Clearly, Eqs. l lS21bt and IS21fl are identical. Furthermore, Eqs. l lS21cl and IS21el contain 
the same information, as can be shown by using Eqs. IS16t l lS14bt . JS17I . and IS21at . The remaining constraints yield the 
equation. 


+ ahoaMO + (“mo + (A:^ - 

+ 0^m£’d{0^jVfO“Mo(“MO “mO^ “mo“mO (“mO^ — Q C^)} = 0. 


(S22) 
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By inserting Eqs. )S13> , IS16t into Eq. )S22t , we obtain a nonlinear equation in O for every value of k. This equation can 
be numerically solved, at least in principle. Q can be then used as input to Eqs. lS14al , )S14bt , IS171 , IS19al , IS19bl , and 
IS21a> - tS21fl to solve for the electromagnetic modes. As we shall see, the most important qualitative feature of the solution 
of this problem is that the magneto-SP fields acquire tangential components (see Eqs. <S14al and IS21bH which are absent 
when g = 0 [I3,[3l, that is, in the absence of an external magnetic field. Eq. )S22t is identical to Eq. (45) in Ei upon carrying 
out the substitutions exx = Czz = £d- 


E. Perturbation expansion on g 

Chiu and Quinn reported a dispersion relation Cl vsk hy numerically solving Eq. IS22> for a very similar setup to the one 
of our interest. However, a detailed description of the resulting electromagnetic modes was not presented in that work. As 
explained, one may in principle solve for the profile of the electromagnetic modes once this dispersion is known. However, 
a numerical attempt at the problem using a standard nonlinear solver yielded spurious results for the modes. 

Since the solution to the SP problem with no magnetic field (g = 0) is a weU-known textbook result, and g is anyway 
much smaller than in a realistic setup, we may use a perturbation expansion of the equation in powers of g. Our goal 
is to obtain the electric fields up to 0(g), so that we can compute the magnitude of the exciton-SP coupling to that same 
order. To accomplish such objective, we first need to solve for Q as well as the coefficients up to 0(g). As we shall see, 
however, knowledge of requires information about Q up to O(g^). Once this is done, we simply Taylor expand the fields 
ISlSat , 9S18bt , )S20at , )S20bl and collect the results according to Eqs. <S19al -l [S20bH . In retrospect, the original problem 
we faced by trying to directly solve Eq. IS22> originated from the fact that we were using the very small O(g^) corrections to 
Cl as an input to solve for the 0(g] electromagnetic modes. This requires an accurate solution of Cl, which is complicated 
by the highly nonlinear dependence of Eq. jS22l on Cl. As a future consideration, it might be worth exploring numerical 
methodologies to attack this problem beyond the perturbative regime, although for our purposes, the latter suffices. 

Even though the algebra below seems involved, it is straightforward to derive using a symbolic algebra package such as 
Wolfram Mathematica(c). 


I. Solving for Cl 

Given that the right hand side of Eq. )S22t is zero, the polynomials at each power of g must each vanishing identically. 
To start with and as a consistency check, at zeroth order in g, Eq. )S22t becomes 


[2cr^0 + 2oido£^mo] ^ i^Oido) + ttmO ~ 

[2a^o + 2a^oO;mo](a:dofmo + a:moedo) = 0, (S23) 

where the 0-subscripted variables denote the corresponding functions in Eq. fS6j . IS13> , jS16> taking Cl - Oq, 

of 

£m 0 — £oo ~ } (S24a) 

ado = (S24b) 

«mo = \J k’^-CmoCll, (S24c) 

where Clp = ^. Eq. jS23l implies that either 

^^do ^cx^QOCjfiQ — 0, (S25) 

or 

ado£mO + OCmO£d - 0- (S26) 

The condition in Eq. IS25I requires that a^o = -o^mOi contradicting the very nature of the SP solution we are looking for, 
where a^Q, amo > 0 represent evanescent fields. However, the condition in Eq. )S26t is simply the standard equation for the 
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dispersion relation of a SP at the interface of an unmagnetized MO sample and a metal film [3, [Hi- It can be readily solved 
yielding, 


k — Oq 


£mO£d 
£m0 + 


(S27) 


or more explicitly, 


Qo - 


1 


[ea + e^)- ^[edO,l +k'^[ea + 




(S28) 


At short k, the (linear) dispersion is very light-like, Oq = and at large k, it plateaus to Oq , corresponding to 

y^d V^<i+£’(X) 

collective charge oscillations in the metal (see Fig.[S2). 

Moving on to the terms of Eq. I IS22I yields an equation of the form = 0 where /{ado, ocmo) is a nonzero polynomial 
in ado and amo, implying that = 0. This result can be quickly derived as follows: the g^ terms stem only from the power 
expansion of a'^Q, a'^Q, am, and Cm- Some g'^ contributions from are proportional to Di but some are not. 

Regardless, the ones from come with the opposite sign to the ones from a'^^, hence, they vanish identically as Eq. 
IS22I is symmetric in and On the other hand, every g^ term for am and Cm is strictly proportional to Oi requiring 
Oi = 0 for the all the g^ terms to cancel. Hence, the lowest order correction to Oq of Q. arises at O(g^), that is. 


O =: Oo + g^02. 


(S29) 


As mentioned, we are solely interested in the calculation of the electric fields in each layer up to 0(g). However, as we shall 
see in the next subsection, these corrections depend on 02- To obtain this coefficient, we expand '^P 

to O(g^), but not beyond that. 


MO ~ “rfO ± SO^dl + g^{OCd20 + 0Cd22^2), 

(S30a) 

<^m ~ “mO + g^ECmZzkiz, 

(S30b) 

Em ~ EmO + g^EmZZ^iz, 

(S30c) 


where the coefficients in the expansions take the form. 


“d 20 

“d 22 

0^m2Z 

£m22 


iklo 

2 v/^ 


kipEd 

«do 


klpEj 

^mO 


202 



(S31a) 

(S3 lb) 
(S31c) 
(S31d) 
(S31e) 


Notice that in order to ultimately solve for 0.2, we have separated the g2 terms into two categories: those which are pro¬ 
portional to O 2 (g^adzzkiz, g^ocmzzkiz, g^Emzzklz) and those that are not {g^ad 2 o)- We substitute these expressions into Eq. 
IS221 and collect the g2 terms, which ought to cancel. The manipulations yield in a linear equation for Q 2 which gives. 


^o^mo(3ed-emo) 


02 = 


^EdiEd - Emo) [o|(2o:2^em0 + ) + n^edEmoiEd - Emo)] 


where we have also used Eq. jS26l in the form of amo-- to simplify the final expression. 


(S32) 
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Figure S2. SP dispersion energy as a function ofk = |fc| for the two-layer (metal-MO dielectric) setup, assuming g = 0. Plots generated using 
Eq. IS28t with Drude parameters for Ag, Coo ~ 4, (i)p ~ 9 eV, and varying the dielectric permittivity e^. Since O = Hq + O(g^), the plots are 
correct up to 0(g), which Is our perturbation order of interest. Notice that the dispersion curves start-off linearly (llght-llke excitations), 
but plateau to constant values at large wavevectors (charge oscillations in the metal), which become smaller as the dielectric permittivity 
increases. 


2. Solving for t^^ 

Next, we expand the coefficient which denotes the contribution of the -i- fields in Eqs. IS19at and IS19bt (t^Q can be 
subsequently found using the constraint in Eq. JS21al ). 


ttio- tdo + gtdi- (S33) 

Here, tdo and are unknown, but we shall solve for them using the boundary condition on the perpendicular electric fields 
(see Eq. jS21c» . We expand and Ezm in Eqs. jS14b> and IS171 up to 0(g^)El 


where. 


^z,MO ~ ^z,dO i gEz.dl + 8 Ez^dZr 
Ez,m ~ Ez,mO + 8 Ez,m2r 


(S34a) 

(S34b) 


Ez,dO — 
Ez,mO — 
Ez,dl — 
Ez,d2 — 


ik 

“do 

ik 

ttmO 

fcQo 

i{adoad2ok+ adoOtd22kkl2-^2kD.on2ed] 


^z,m2 ■ 


ikQoQ2£i 


a 


3 

mO 


(S35a) 

(S35b) 

(S35c) 

(S35d) 

(S35e) 


We plug these expressions into Eq. IS21cl . At zeroth-order in g, we get CmoEzmo = ^dEzdOt which is equivalent to Eq. IS261 
and does not give us information on tdo or fdi (this indeterminacy ultimately reveals why we need expansions up to 0[g^) 


1 Here, as opposed to Eqs. [pOi]-(llo 3 we have not separated the determined in Eq. [H. 

terms in and Ezm S^EzmZ) into those contributions 

which are proportional to 02 and those which are not, as O2 has already 
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to obtain fields up to 0(g)). At 0(g) we get the intuitive result, 


tdo 2 ’ 


(S36) 


and at O(g^) we obtain, 


^dl - 


^z,d2^d ^z,m0^2^m22 + Ez,m2^m0 

^^zdl^d 

, ( 1 60:^0+ Sal^n^gOzepoemo + a^oOp + 

=rrfi 


(S37) 


which, together with Q 2 in Eq. IS32t , can be readily evaluated with zeroth-order parameters. We notice that tdi is purely 
imaginary-valued, so we have written it in terms of a purely real-valued t^i. Given Eq. I IS33I , it is clear from Eq. lS21al that 


^MO ~ girfi- 


(S38) 


3. Collecting the expressions for the fields 

We now have all the ingredients to evaluate Emo and Bmo (see Eqs. IS19a> - tS19bt ) up to 0(g), 


EMO-Edo + gEdi, (S39a) 

Bmo ~ Bdo + gSai, (S39h) 

where, at zeroth-order, in the absence of external magnetic field, we have the standard SP mode which is a transverse 
magnetic (TM) mode fl3.[l5|] , 


Edo = (l,0,—)r7oe"“‘''’", 
i adQ> 



icrfQo 

^do^ 


(0,1,0) 77pe"“‘«>^. 


=Bdo 


(S40a) 


(S40h) 


In going from the first to the second line of Eq. jS40ht . we have used Eq. l lS24ht . Eqs. IS40at - tS40bt feature an elliptically- 
polarized electric field with no tangential component, and a purely tangential and imaginary-valued magnetic induction. 
The opposite is true for the first-order correction: it consists of a purely tangential and imaginary-valued electric field (recall 
that fdi is purely imaginary, see Eq. 1S37I ) and an elliptically polarized electric field with no tangential part. 


Edi - 


/no(4Trfiv^-Qpz) 

2 UdO 


( 0 , 1 , 0 ) 770 


Bdi - 


1 


=Edl 

^ocdoT^di VEd + Go(1 - arfo z),0, ik(4Tdi - noz)| r]oe 

=Bdl 


(S41a) 


(S41h) 


In deriving Eqs. IS39a> - IS41a> . we have used the 0(g) Taylor expansion for the fields (see Eqs. IS14at - lS14bt . iSlSbt l. 
Notice that besides the exponentially decreasing dependence of the fields, we also obtain a polynomial contribution in z. 
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We already computed some of the relevant expansion coefficients in Eqs. IS34al -l [S35cl : the remaining ones that we used 
are, 


T7+ 

^8,MO ~ ’ 

“do 

(S42a) 

n-t _ * ^0 

(S42b) 

. iElo^d 

a,oC ^2a2„c’ 

(S42c) 

± ^ 

^z,MO ~ ~ rr , r’ 

“dOd 

(S42d) 


as well as those for the plane wave components, 




, gi(fcr-enof) + 

=m 


iClQ 


(S43) 


Finally, given Emo, Em can be readily obtained from the boundary conditions for the fields at z = 0 (see Eqs. jS21al - iS21fH 
as well as the original ansatz for their functional forms (see Eq. jS20al - fS20bH , 


Em ~ EmO + S^mlt 

where, at zeroth-order we have the standard SP electric field as if there were no magnetic field present, 


(S44) 


= )77oe“"■“^ 

—^mO 

Brno = ( 0 , 1 ,0)J70( 0 , 1 , 0 ) 7706“'"“^ 

OCmOC OCiIqC 

'- , -' 

=BmO 

Here, we have used Eq. jS26l in both lines. Similarly, the first order correction to Emo is, 

Emi = (0,1,0) 770 

“do 

'- , -' 

=Bml 

Bmi = ^^[oo-a:do\/^.0,7A;^^]77oe“'"“^. 
“docl I 

=Bml 


(S45a) 


(S45b) 


(S46) 


F. Quantization and normalization of modes 

The previous section shows how to (perturbatively) compute the frequency Q, the electric field, and the magnetic in¬ 
duction of an SP mode with wavevector k. Note that, so far, we have invoked an arbitrary normalization (setting 1?^ = 1). 
To fix this, we first ought to compute the energy associated with the unnormalized modes. Consider placine electric field 
amplitude .e/fc into the radial component of the k mode. The energy in this mode is quadratic in the field^ lla . [l^ . 


^ Note that even if the electric field is much larger than the magnetic in- tributions to the energy density are on the same order of magnitude 

duction, the prefactors of cq and weigh them in a way that their con- 
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Hsp.k--Z / eoZ 


d(i) 


{ErAE]i +—mif 

■' MoM 




(S47) 


where i,je{r,6,z}, and electric field and magnetic inductions in each k mode (throughout z) are conveniently written as, 


E - Q{-z)Em + Q[z)Emo> 
B = Q[-z)Bm + &[z)Bmo^ 


(S48a) 

(S48h) 


and {E)i, {B)i denote the i-th components of the respective fields (which include the plane wave exponential factors, see 
(SUD-illOb)). InEq. tS47t . we have absorbed the electric field units into sik, so Ei is taken to be a unitless quantity. The 
integration f dV is carried out over all 3D space. 

Assuming a finite size box of in-plane area S and infinite perpendicular dimension and plugging in Eqs. IS19at , IS19bt , 
IS20at , and l lS20bt into Eq. IS47t , we obtain. 


/ dzUoE — :r — iE]*(E]i + —mi\^] 
2 / 7-00 ^ i d(u J poft > 


d{(jDe*j{(jD)] 


S\.S^kf V" ■sp , r Vi-* fpl' 'l*P^ , ^ rnf )*dS 

2 '-‘•mo’ ‘mo ‘■O Z^‘-MO,ij'-^j,MO’ i.MO’ ^i.MO 


i 7,5e{+,-l 


dze ■’■“moD 




E 


d{wem{M)) 2 1 ID |5 

£0 'Z \Eim\ “I \Bim\ 

d(t) fiQfi 




- S\^k\^J^ E ^^MO^ ^MO ^oE^MO.O'^^J.MO^ ^^MO^ ^ImO 

i jm+-] ‘ i Bob 


2{aj + as) 


+S\s^k\^Y. 


diOJEmm 2 , Ip ,2 

^0 , \^i,m\ 1-0/,ml 

do) fiQiJ 


1 




where, following [ii] , we have defined the vertical mode length as. 


(S49) 


i-fc-E E ^MO ^‘’Yj^MO,ij^^],MO^ ^f,MO^ ..^^1,M0^*^UMO 

i 7,5e{+,-} j BOB 


eoia* + as) 


+E 


fo- 


d{( D£m[OJ)) 
d(i) 


+ \Bim\ 
BOB 


1 


eoiam + am) 


(S50) 


Since Ei is taken to be unidess, Lfc effectively has units of length. Physically, Lk defines a mode volume SLk with a quantized 
amount of energy corresponding to the frequency Wk- Importantly Eq. <S49I is quadratic in the electric field amplitude s^ki 
which implies that each k mode corresponds to a harmonic oscillator. If we wish to quantize the energy in quanta ofoJk, 


U *, * A 

^^SP.fe = -t aj.afc), 


(S51) 


we can define at so that. 


- 



(S52) 
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Other multiplicative phase factors in this amplitude definition (U(l) gauge choice) do not affect the quantization. Promoting 
the complex amplitudes to operators, a/t and a*^ —► with [ajt, = 1 , 


(i){k) 4- f 

Hsp,k - + a^aic) 




(S53) 


To summarize, we have normalized each k mode in Eqs. )S48a> and jS48ht hy associating energy quanta Wfe to a SP excitation 
in such mode. Finally, the final electric field and magnetic induction are superpositions of amplitudes in such modes. 


S^Y.s^kE{k]. 

k 

k 

Promoting these amplitudes to operators (in the Heisenberg picture) and using Eq. jS52l . 


(S54a) 

(S54b) 


?^V2eoSLfe 


u)[k) 


auE{k] 


„ , (i){k] 


(i»[k) 


SS(.r,t) = y ,, 

V V 2eo SLfe 


(W(fc) 


aicBik) + 


2eo SLk 
a)[k) 


alE*{k), 

(S55a) 

a\B*{k), 

(S55b) 


2 eoSLfc 

where we have used the fact that S and SS are real valued to write the operators in a more conventionally symmetric form. 


G. Perturbation expansion of L ^ 

The fields deduced in Sec. [T^are only correct up to 0(g), so it is important that we keep Lfc up to that order too. Impor¬ 
tantly, since (*^ mo) re - ig and Ega oc g, the lowest order contribution to Ljt appears at O(g^). Hence, we do not need to 
take the anisotropic effects of cmo into account, and we can set g = 0 in Eq. {ST), i.e. ^mo ~ where D is the 3 x 3 identity 
matrix. Therefore, the calculation for Lfc L^o reduces to that of the standard SP mode in the absence of external magnetic 
field (a^o and amo are purely real) El], 


Lko = y 


^O^dlEidol + l^cdol 

pop 


1 


26’ocr^o 


eo- 


d(wCmo(^d)) 


dco 




'■do 

~£mO 1 
■ + ■ 


2ado 


dicjemoio))) 


ndo 2oCffiQ 


d[(Demo{(»]) 


d(o 


. ^0 


d(!) 

PmO ~ £d 
£m 0 


l^cmol + |fi(,mol 

£o pop 


2 eoamo J 


I ado I 


2am0 


■^mO 


(S56) 


In this derivation, we have used p = 1, eopo = c Eqs. IS27> . IS35a> . and lS35bt as well as the expressions for (S24b), and 
(S243M 


The final result in Eq. [S56I is twice of what is reported in the Supplemen¬ 
tary Material of Ei . We believe our derivation has the correct prefactors; 


however, the use of either result gives the same order of magnitude of the 
effects we are interested in. 
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z 



Figure S3. Three-layer (metal-MO dielectric-organic) setup. We are interested in the (magneto)-SP modes arising at the dielectric-organic 
interface (z = a) upon application of a perpendicular external magnetic field. 


II. MAGNETO-SPS IN A THREE-LAYER SETUP 


We shall now adapt the results from the previous sections to the situation where the MO layer has a finite height a, and an 
organic layer of isotropic dielectric Corg is placed on top of it (see Fig. [^. As far as we are aware, the resulting expressions 
for the corresponding magneto-SPs have not appeared before in the literature. 


A. Electromagnetic modes for each layer 

1. Organic layer (z> a) 

Just as we did for the MO and the metal layers (see Eqs. ISIOI and IS15U . Eq. for the organic layer can be expressed in 
matrix form Mlorgitorg = 0, 


£m ^^z,org ® ^ ^ kkz.org 


^k,org 


0 

0 em-n~^{k^ + klg^g) 0 


^9,org 

= 

0 

kkz.orgEl~^ 0 Cm-Er^k^ 


_^z,org . 


0 


(S57) 


Importantly, we do not fix the radial component Ek.org to 1 because of the boundary conditions at the organic crystal and 
MO interface at z = 4E The corresponding secular equation for the decaying field for z> a, with kz.org - iocorg yields 


^org — 


\Jk'^ - O^e 


org, 


(S58) 


In our (arbitrary) normalization before quantization, we may set only 
one of the field components in one of the layers to 1, and our convention 
is to choose Er,m = 1 as in the two-layer case. The rest of the fields are 


not arbitrary and satisfy the wave equation Eq. fS9] as well as boundary 
conditions at each of the interfaces. 
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as expected (see Eq. IS16H . Finally, just like in the metal layer, Eq. IS571 tells us that 


^z,org — Er,orgi (S59) 

aorg 

hut does not inform us about the tangential component Eg^oj-g (nor about Ek^org)i The missing components will be deduced 
by matching the boundary at z = a. Altogether, the fields in the organic layer read like those in Eqs. ISlSal and ISlBbt , 

Eorg — EgrgTje 

— iEif^Qfg,Eg ffi,Ez,m)TI^ I (S60a) 

Eorg — Eorgl]^ 

— —^{(XorgEQ^orgj~^^Ez,org~^orgEk,orgj^^Eg^Qfg)Tj6 , (S60b) 

denoting exponentially decreasing fields. 


2. MO layer fa> z>0J 

For the MO layer, we translate Eqs. ISlSal and ISlSbl to this setup, 


= (l.-Efl A/fni , 


^MOl ~ 


e.MOl’ z.MOl-’ 
a-z 


-I 

t) 




7± 
"z.MOl 




(S61a) 


(S61b) 


where and indicate exponentially decreasing fields {kz = and by slightly adapting these expressions. 


Emo\ - Emo\E^ 

I7± , 


“ ^^’Eg.Mnt’Ez.MO\^''l^ 


^MO] ~ EmO\E^ 


-I 




(S62a) 


(S62b) 


where and denote exponentially increasing fields {kz = In previous sections where MO was consid¬ 

ered to fill up aU the space z > 0, the latter fields were not considered, the reason being that e“Mo^ was unbounded as z ^ oo; 
this is not the case when the largest value of z is £20. Hence, the analogous expressions to Eqs. IS19al and IS19bl are. 


Emo - Emo[ + Emo\’ 

BmO - Ij ^MOl h ^MOl ^MO\ ^MOf 


(S63a) 

(S63b) 






e,Moi 

+ 


7± 

"e,MO\ 


identities are easy to check as is associated with kz = and E 


were derived for kz - ioL%[Q- 


"MO\ 


'e,MO ^z,MO[ ~ ^z,MO\ ~ ^z,MO- These 
with kz = -ia^Q, but Eqs. IS14al and IS14bl 


5 A more intuitive way to describe this situation is that and ggijs starting from z = a going downwards, 

denote fields that exponentially decrease starting from z = 0 going up¬ 
wards; similarly, and B^jq^ describe exponentially decreasing 
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3. Metal layer (z < 0) 

Finally, for the metal layer, all the expressions we derived in Sec. Uhold, in particular Eqs. jS20al and IS20bt : we still 
assume Ek,m- 1- 


B. Matching the modes at the boundaries (z = 0 and z = a) 

Given this prelude, the analogous boundary conditions to Eqs. iS21al -l lS21fl for z-0 are. 


1 = t 


MOl ^MOl '•MOt ‘MOt’ 

Ee,m - ^MOl^e.MOl + ^MOp^e.MOt + 


+ h. 


• + t. 


CmEz.m - ^di^MOl^z.MOl + ^MOl^z.MOl + ^MO]^z,MO] + ^MO]^z,MO]^’ 


^E^ fyi + itXifi — lc{t 


+ p+ 
MOl^z.MOl 

MO'-'-MOl ~ ‘MOt 


+ ^MOl^z.MOl + ^ 


+ tj. 


]- ia 


MOt^zMOt ' '■MOf^z.MOt 

), 


MO^^MOl ^MOt 


^B.m - ^mOI^BMOI ^MO[^B,MO[ ^ ^MO\^6,MO\ ^MO\^B,MO\’ 


whereas for z - a they are. 


(S64a) 

(S64b) 

(S64c) 

(S64d) 

(S64e) 

(S64f) 


_++--+! - 1 
Ek,org - ^MOi^ + ^MO\^ + ^MO\ 7+ ^MO\ 7^’ 

A A 

^B.org - ^MOl^e.MOl^^ + ^MOl^B.MOlX + ^mO]^ 6 ,MO] + ^MO\^B,MO\~’ 

A A 

orgEz,org - + ^MOl^z.MOl^ + ^MO]^t,MO] + ^MO]^z,MO] 7^ 


OtorgEg^org- | “mO^MO1-®0,MOI^ ^ ^M0^M0[^B,M0['^ 

+ + + 1 _ _ _ 1 ^ 

+ 1 O^MO^MO\^BMO\~j^ ^ ^MO^MO\^B,MO\~j^j’ 

kEz,org - iaorgEr,org - + ^MOl^z.MOl^ + ^MO]^tMO] + ^MO]^z,MO] 7^ 

A A 






1 


‘■MO 




^B.org - ^MOl^g,MOl^ ^MO[^e,MO\^ ^MO\^B,MO\ + + ^MO\^B,MO\ y- ' 

A A 


(S65a) 

(S65b) 

(S65c) 

(S65d) 

(S65e) 

(S65f) 


Here, j- = embodies the vertical thickness dependence of the problem. Eqs. IS64a> - tS65fl can be manipulated to 

yield an entirely analogous expression to Eq. <5221 . This resulting expression can then be perturbatively expanded in g and 
the analogous procedure of Sec. ITElfollows. with the caveat that the algebra becomes much more laborious. We bypass the 
latter by automatizing such work with Wolfram Mathematica(c). We summarize the results in the following subsections. 


C. Perturhation expansion on g 

1. Solving for and 

It is clear that provided that the frequencies Oq and Q 2 are modified accordingly (one finds Oi = 0 again), the expansions 
for the dielectric constants and the wavevectors can all be recycled from Sec. [T^ On top of these quantities, following Sec. 
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Figure S4. SP dispersion energy as a function ofk= |fc| for the three-layer (metal-MO dielectric-organic) setup, assuming g = 0. Plots 
generated using Eq. IS69I with Drude parameters for Ag, Coo ~ i, ti>p ~ 9 eV, and varying the dielectric permittivity as well as that of 
the organic layer. Since O = Oq + O(g^), the plots are correct up to 0(g), which Is our perturbation order of interest. Just like with Fig. 
IS2I , notice that the dispersion curves start-off linearly (light-like excitations), but plateau to constant values at large wavevectors (charge 
oscillations in the metal), which become smaller as the dielectric permittivity increases. 


ITEI we ought to expand ctorg and x~ up to 0[g^) in Eq. <S58I . 


ttorg ~ ttorgO + S ^org22^2> 


ttorgO — ^gSorg: 

^O^org 

ttorg22 — > 

ttorgO 


(S66a) 

(S66b) 

(S66c) 


X~~Xoil+agadi-g 




(S67) 


but we only need t^^^ and up to 0(g), 


^M 0 l~ ^dio±gtdii, (S68a) 

^MOt ~ 

Collecting the zeroth-order in g contributions in Eqs. IS64a> - tS65fl , we derive an implicit equation for Oq which coincides 
with the standard textbook result for a three-layer system in the absence of an external magnetic field Q, 


2 _ { ^dO^mO + ttmOCrf ) I O^dO^org + OtorgO^d ) 

I ~ / V (^d^^org ~ <^org0^rf / 

where a^Oi oLma, ocorgo, and Cmo are all functions of Qq- This equation is the analogue of Eqs. )S27t and )S28I for three layers. 
It is not possible to explicitly solve for Oq- but one can readily compute it numerically from such implicit equation (see Fig. 
[S4). The formula for Q 2 is too long to display it here and is, anyway, not relevant on its own. However, we make use of it to 
obtain our final results. 
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Next, we compute the coefficients in Eqs. )S68a) and IS68bt in analogy to the two-layer case in Eqs. <S36I and jS37l . At 
zeroth-orderin g, 


tdlO - 


td\0 - 


^mO^d ^dO^mO 


(S70a) 

(S70b) 


while their first-order in g corrections are, 


tdii = ocdi[aa^M^xl-'^'>£mo 

~^d0^^0 (^mO + ^orgo) + ^^mO^d + ^mol + ^^mO^d ~ ^^mO^mO + ^^orgO^mO ~ ^mol 

-Harfol^a^oUo + + OCmoiaaorgO + 1) [xl(-^d + ^mo) - + e'mo] + 2a orgO^^mo] 

~^mO^Xo~ ^^^d(<^^orgO + l)j / 

{4amo£’d[lo(amo - ado)iocdo - aorgo) + {ado + amo){ado + aorgo)]] (S71a) 


and 


td\l — ~ aorgO [a{cCdo + t^mo) "t 2] -I- CKXdoamO + ado ~ amo] 

~amo^d {ado + afno)\a.{OLdO ~ aorgo) ~ 1 ]] 

~ad\ {ado + ®mo) lai<XdO + ®orgo) + 1] {adO^mO + j ! 

{4amoerf[lo(“mo - ado){ado - aorgo) + {ado + amo){ado + aorgo)]}- (S71b) 


Hence, at zeroth order, -i- = 2tdio and + t^Qj = 2td]0y and these total coefficients for exponentially decreas¬ 

ing and increasing fields become identical to the textbook results for SP modes in the three-layer setup in the absence of an 
external magnetic field. 


2. Collecting the expressions for the fields 
For reference, the zeroth-order fields are given by. 


EorgO — {Ek,orgO>^>^z,orgo)tJoS > (S72a) 

Bor go = (0, Be,orgo, 0)770 (S72b) 

Edo = 2 fdjo (1,0, E,_dio)Vo 6"“''''^ + 2 tdf 0 (1,0, Bz,rf^o)??oe“'«>^ (S72c) 

Bdo - 2frfjo(0,Be,rfjo,0)7706“““''’^ -i-2tdto(O,B0,dfo.0)7706““'“^, (S72d) 

B^o = (l,0,£^,mo)?7oe"“”'“^ (S72e) 

Brno = (0,Be,mo,0)poe"“"-'>^ (S72f) 








25 


with each of the components being, 


ttmolXo + - “4o(lo “ IJ^mO 

Bk.orgO — „ ) 

2am0l0fd 

(S73a) 

^ _ ik[ado{xl + '^)Cmo + oCmo(£d-xled)] 

2adoamoXo£org 

(S73b) 

i{2adoamoaorgoXo£org + k^ladoixl + De’mo + amo(£d - xl^d)]) 

Ee,org0 — no ’ 

^‘OCdoOL.mOXO^O t^^org 

(S73c) 

ik 

Bz,d\G — ~Bz,d\G — > 

CtrfO 

(S73d) 

i(k^-a^,) 

Be.dia - Bg^d\Q - „ 

adQ^lQC 

(S73e) 

ik 

^z,mQ — - f 

^mO 

(S73f) 

-D0,mO - ^ 

(S73g) 


The expressions for the 0(g) fields at each layer are also cumbersome; we only show the electric field in the organic layer, 
as it is the one associated with the coupling with excitons (see Eq. )S39aH , 

^org ~ EorgO + gEorgl- (S74) 

The first order in g contribution is 

= (S75) 

where Eorgi is purely tangential and purely-imaginary valued, 

Eorgl ‘0 — ~ ~ ^orgo) + ^orgol ~ Xo~^ 

~^^mo{XQ^^^doi^orgO ~ ^do) + ^orgol + ^(^doi^dO + ^orgo) ~ tt^orgolj/ 

2(Xm0^dlXo^^TnO ~ ^do)(.^dO ~ ^orgo) + i^dO + ^mo)i^dO + ^orgo)]- (S76) 

Compared with Eq. IS41at . the z-dependence of Eorgi is purely exponential, as the lowest order correction of aorg to UorgO 
is 0 (g 2 ). 


D. Quantization and normalization of modes 


Equipped with these results, we use Eqs. IS64al - R65fl to compile expressions for the fields in each layer. First, we aim to 
compute the vertical normalization length Lfc for each mode. In analogy to Eq. IS5QI , we obtain, 




^0^org\Ei^org\ “t" lEi^orgl 


^0 J a 


+ E E Y ^MOv 

ij r,(5e{+,-} li,!/£{(,11 

1 


1 


eo(£MO,ijE^-^j^o) Eij^q+ 

MOM 


ra 


£o Jo 


■E 


eo- 


du) 


\Bim\ “I \Bim\ 

Mop 


^0 J —CO 


(S77) 
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where sgn(|) = -sgn(i) = 1 denote the exponentially increasing or decreasing fields, respectively. At 0(g), Lfc ~ L^o as in 
Eq. 1S56I , and the anisotropy of the MO layer is unimportant. Carrying out the integrations explicitly (and taking care of the 
possibility that a do is complex-valued), 


Lko - Y. 


^O^org\^i,orgo\ l^z,orgol 


2eQ^OL.orgO 


+4^ tdtol 

i 

+4^ Itdiol 


£o£d\Ei,d]or +- \Bi4\or 

MojU 

£o£d\Ei,diof' -t- IBi^diof' 

^^o^^ 


g2!Ra,joa _ 

leo^ado 

\ — g-29tarfo® 

2eo5Rado 
1 


+ ■^4V(rrffo)*(fdio) ^oY^'iiEi^d\o)*Ei4\o +- [Bi4^o)*Bi4io 

i] L / Bo^J■ 


a 

— I -i-c.c. 
eo 


-E 


eo- 


d{(D£m(.(i))) , 


1 


, '\Eimo\ + \Bimo\ 
act) Hold 


1 


2eoRamo ’ 


(S78) 


One may numerically compute Lt by plugging Eqs. IS70a> . lS70bt . and IS72at - (S73g) into Eq. I IS78I . We do not display the 
resulting analytical expression, which anyhow, is lengthy and not particularly illuminating. 


III. EXCITON-EXCITON AND EXCITON-SP COUPLINGS 


In Secs. Uandini we solved for the electromagnetic profile of the magneto-SP modes in a two- and three-layer setup. 
We are now ready to describe the organic superlattice. We regard the latter to be either “embedded” in the MO layer (in 
the two layer setup) or in its separate third layer (in the three-layer one). As explained in the following paragraphs, the 
superlattice consists of a monoclinic array of organic aggregate nanopillars. For simplicity, we take each of the nanopillars 
to be a rectangular parallelepiped of volume WxWyWz (here, W; is the width of the nanopillar along the i-th axis). If the 
nanopillar density is pnp, it contains Nnp - PnpWxWyWz chromophores. Furthermore, the three-dimensional positions of 
the individual chromophores constituting each nanopillar are denoted. 


rms= {mx8xX+my8yy)+ s^z, 
= rm 


where d/ is the spacing between chromophores along the i-th direction. 


(S79) 


A. Dipolar couplings between nanopillars 

We shall first study the energetic contribution due to the excitons alone. We model this as, 


Hexc = Y‘d^n(yh(dn+ E (/nn'^^Vn'+ h-C.), (S80) 

n^n' 

where is the bare energy of the n-th collective nanopillar dipole, which we take in the ideal case to be = (h. The exciton 

hopping amplitude between the n-th and n'-th nanopillars is approximated as a near-held dipolar coupling. 


Jnn' = 


e\rn-r„'\^ 






(S81) 


q 9 1‘n~f ' 

where q = 0.625meV(nm /D ), e„„' - -r",\ ’ where r„ is the average in-plane location of the n-th nanopillar, and we take 

e =: 1, the dielectric permittivity in the medium surrounding the nanopillars. Eq. <S8U implicitly relies on a separation of en¬ 
ergy scales, namely, that the coupling between chromophores is much stronger within a single nanopillar than between dif¬ 
ferent ones. Hence, we start with the collective superradiant nanopillar transitions which scale as = Hms Pms ~ \/^npPn’ 
where is the transition dipole moment of the ms-th chromophore in the nanopillar, the sum is over all ms values as¬ 
sociated with the w-th nanopillar, and - p„, that is, we take the dipole to be equal for aU chromophores within the 
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corresponding nanopUlar. This approximation should provide a semiquantitative description of the dispersion of the or¬ 
ganic superlattice alone. A more refined description would rely on the coupled-dipole method [T^ - [l9|] , but is beyond the 
scope of this work, as this simplified model illustrates the essence of the problem. 

We shall now consider a general two-dimensional monoclinic superlattice with unit cell defined by vectors OD and OC 
depicted in Fig. [S^ For convenience, we temporarily adopt the 0-rotated coordinate system x'y' which, with respect to the 
original xy system, is defined by, 



COS0 

sin0 


-sin0 

COS0 



(S82) 


We win later explain how to obtain 6. We take two sides of the parallelogram {AB and CD) to be parallel to x'. For simplicity, 
the nanopiUars are taken to be rectangular parallelepipeds. Their transition dipoles n„-fi are fixed in the x'y' plane and 
make an angle a' with respect to x' {or: a = a' + 8 with respect to x). Notice that all sites are equivalent. We only account for 
nearest-neighbor (NN) and next-nearest-neighbor interactions (NNN). We classify the interactions as horizontal NNN {AB, 
CD), vertical NNN {AD, BC), diagonal type A NN {OA, OC), and diagonal type B NN {OB,OD), respectively, 


n (1 - 3cos^a') 

-73-’ 

(S83a) 

P [l-3cos^(a'-d)] 

Jv-VB 3 

(S83b) 

2 [l-3cos^(a'-7)] 

diagA IB r2i2^+2Ah^vCOSP+Al\3l2’ 

(S83c) 

^ 2 [l-3cos^(a'-i-(5)] 

JdiagB - VB jA2_2A,A„COSyS+A2 j3/2 ' 

(S83d) 


Here, |^| = /i and the side lengths are AB - CD - and BC = AD - Ap. We have also conveniently introduced the angles 
fi = ZBCD - ZDAB, as well as the following. 


7 = atan 
6 = atan 


Aysin;6 
Ah + AyCOSjS’ 
Ai,sin;S 
Ah - Ai,cos/3 


(S84a) 

(S84b) 


Assuming that all site energies are equal, = o), we may rewrite Eq. jS801 in k space, Hexc = T.kHexc,k> where Hgxc.k - 
(i>exc,k(^\.(^k< where k - kx’x' + kyy' and 


(^exc,k + 2JhCOS{kx'Ah) + 2/i,cos 


kxiApCosp+ ky ApSinp 




A/,-t Aycos^d Ai,sin;6 


+ k 


■y- 


A;j-AyCosj8 Ai,sin;6 
kx' - kyl - 

■^9 / 9 


(S85) 


is the resulting dispersion relation for the excitons alone. 

As explained in the main text, we would ideally like to design a plexciton dispersion which features a global gap in the bulk. 
This requires a superlattice with an “H-aggregate” dispersion along all wavevector directions. Mathematically, this means 
that the dispersion (i)gxc,k should be a maximum at fc = 0. It turns out that this is not possible in a rectangular lattice {fi- j), 
as the resulting J-couplings {Ji < 0) arising from the geometric constraints end up dominating the H-couplings {Ji > 0) at 
least along one direction, yielding a minimum, or at best a saddle point for ojgxc.k at fc = 0. Hence, we proceed in a more 
systematic fashion. Taylor expanding Eq. 1S851 up to quadratic order in fc,;, 


























28 



Figure S5. Geometry of two-dimensional monoclinic superlattice of organic nanopillars, x'y' denotes a temporary Cartesian coordinate 
system which is rotated at d from the original xy system. The collective transition dipole moments of the nanopillars make an angle a! 
with respect to x' (or a= a' + 0 with respect to x). The horizontal and vertical distances A/j and Ay, together with the angles f, 7 , and S 
fully define the superlattice. 


<^exc,k ~ I 

where the constant offset is obtained by evaluating otgxc.k 


Mx'x' 

Mx'y' 

^x' 

_ My'x' 

My'yl 

kyf 


at fc = 0, 


(S86) 


^eff ~ + 2J)i + 2/y + 2Jiiiiigj{ + 2JiiiiigB- 


t-fp^p A/r., - I ^^^cxc.k 

nere, 2 dk,,dki, 


denotes a Hessian matrix, which can be readily diagonalized, 
fc=o 


^x'x' 

AJx'y' 


Sxx' 

Syx' 

mx 

0 

Sxx' ^xy' 

My'x' 

Mylyl 


Sxy' 

Syy' _ 

0 

my 

Syx' Syyl 


where S/y' is a unitary matrix, and ultimately yields the result, 


(S87) 


(S88) 


o^exc,k~^'eff+mxkl-kmyk?y (S89a) 

= (heff + ZJxCoskxAx + 2/yCosfcyAy, (S89b) 

where fc; = Six' kx' + Siyi kyi, and for our simulation, we choose (arbitrary) effective unit cell dimensions Aj such that Ji- - ^ 
for i - X, y, and f + 2Jx + 2Jy; this identification renders Eqs. )S89at and )S89bt equal up to quadratic order in ki. 

We have thus arrived at a very convenient expression; Eq. jS89bt shows that the oblique lattice renders the same long- 
wavelength physics as a much simpler rectangular lattice with only NN interactions. This approximation is insightful in that 
it exposes the physical origin of the global gap; it also remains valid for our purposes as the topological phenomena of our 
interest occurs at small fc,-. 

For a fixed value of A^, we Monte Carlo sample through the parameters r = ^ e [0,2], a e [0,7i], and p e [0,^^] and record 

those which yield ntx, my < 0. We observe that only ~8% of the parameter space satisfies the H-aggregate condition we are 

looking for. One such set of parameters is a' - 1.19(~ 68.5degrees), p - 0.23(~ 13.1 degrees), r = 1.2, yielding [mx, my) - 

2 

-(1.49,0.44)/oA?, where Jo = ^ sets the energy scale of the dipolar interactions. The associated eigenvector matrix is, 
n A, 
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^xx' ^yx' 


-0.91 

-0.41 

^xy' ^yy' 


0.41 

-0.91 


COS0 -sin0 
sin0 COS0 


(S90) 


where we obtain 6 - 2.72(~ 155.7°), which defines the angle of rotation of our temporary coordinates x'y' with respect to 
the original ones xy; then, a - a' + 6 - 3.9(~ 223°). Now, each nanopiUar has a collective transition dipole moment value of 
po = \/Nnp\p„\ - \/Nnp X lOD = \/WxWyWzPnp x lOD. Choosing the nanopillars to be separated from one another by A/,, 

TlWx^V^v^V^zPnD 7 ^ 

we get a value for the energy scale of Jo - -p^ x 100 D . Taking p„p - 37 molecules/nm , we obtain. 


A == 3.45 X loVeVx G 

rriy , 

Jy^ - ^^1.02x 10^ meV x (y 

Ay 

where (i - are dimensionless ratios which govern the effective dispersion of the superlattice. By choosing the phys- 

ically reasonable parameters A/; = 100nm, Wx = 7.5nm, Wy - 50nm, - 70nm, and = 2.15 eV as well as the effective 
simulation parameters IS.x-l^y- 50nm, we obtain Ay = rA/j = 88nm, (x = (y - 0.105, Jx = 362meV, Jy = 107meV, and 
(Dgyy = 3.09 eV. We emphasize that Jx and Jy can have arbitrary values which depend on our choice of parameters Ax and 
Ay. The latter set the spatial resolution of our real space simulations, and hence, the size of the systems we can computation¬ 
ally study. These simulations are carried out in order to calculate edge states as well as understand the effects of disorder. 

We note that when choosing parameters, we need to make sure that (a) the nanopillars do not yuxtapose each other and 
(b) the number of chromophores in the organic layer is large enough to achieve considerable coupling with the SPs (see Eqs. 
IS961 and 1S99I ). Condition (a) is easily checked computationally and graphically. With respect to condition (b), we note that 
the surface area of the ABCD parallelogram in Fig. [^is A/jAySin/i. It contains two nanopillars of surface area WxWy. The 
surface coverage fraction of the organic layer is hence. 


2WxWy 

/=-(S91) 

A;,AySin/l 

For our chosen parameters, / = 0.38. In general, we need both / and to not be very small (/ >~ 0.2, Wz >~ 40nm). 

In the next subsections, we shall work with the effective rectangular superlattice of Nx x Ny nanopillars (where Ni is the 
number of nanopillars along the I-th direction) instead of the original monoclinic one. This is a good approximation not 
only for the interactions between the various nanopillars, but also for the exciton-SP couplings, as long we use the average 
density of the original monoclinic lattice (see Eqs. I IS96I and IS99) ). 


B. Exciton-SP couplings 

We are now ready to discuss the effective interaction between SPs and a single nanopillar. Consider the dipole operator 
Pms - Pms^^lns + i>ms), where b\,ig[bms) creates (annhilates) an exciton at the ms-th chromophore of some nanopillar. The 
time-independent electric field operator is <?'(r) = Zfc \J 2 e^SLk ^kE'[k) + h.c., which results from transforming Eq. IS55al 

from the Fleisenberg to the Schrodinger picture by removing the dynamical phases (see Eqs. IS18al and lS20al ), i.e., 

E'{k) = Using Eq. jS74t and )S751 , the dipolar coupling between the rath nanopillar and the SP modes is given by 


jL/tw) 

^exc-SP ■ 


m,s 


k,m,s 


(S92) 


where the sum over ms is restricted to the chromophores in the n-th nanopillar. We have also used the rotating-wave ap¬ 
proximation to discard far-off-resonant terms of the form and a^i^crhs- Using Eq. IS55a) , the corresponding coupling 

is given by. 
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c^k,ms — 




(S93) 


where, depending on whether we use the two-layer or three-layer setup results, we make the following substitutions. 



For two-layer setup 

Lk 

Tfco 

Eq. (Ho) 

a 

“do 

Eq. (S24bl 

E(k,Zs) 

EMoik, Zs) = £do(fc) + gEdi (fc, Zs) 

Eqs. lS40al, IS41al 



For three-layer setup 

Lk 

Tfco 

Eq. 1S78) 

a 

^orgO 

Eq. (S66bl 

E(k, Zs) 

^org(^} ^s) — ^orgoi^) + S^orgli^r ^s) 

Eqs. lS72al. IS73al, lS73bl, (S761. 


Eq. <S93I exposes the 3D nature of our problem with its Zs dependence: there is an exponential contribution from 

the SP evanescent field, and even a linear correction in Zj due to ii^i for the two-layer setup. In any case, it will prove 
convenient to derive an effective 2D description for our model. We have two ways to do so. 


J. Mean-field approximation (MFA) 


Using Eqs. IS92I and JS93h 


rrin) _ _ y- 

‘^exc-SP ~ 


(W(fc) 


^ Y 2eoSLj; 


^kPms 


Y,E{k,z Cm n) o( ( 5 ( “ + C.C. 






^ y 2eoSLfc 


(S94a) 


(S94b) 


where we have formally taken z{k) to be an average (fc-dependent) vertical position for the chromophores in the nanopillar 
(we will discuss how to compute this parameter later, see Eq. ISIOIH . made the MFA that ~ y, assumed 

that the dipoles - p„ for all the chromophores in the n-th nanopillar, and defined the collective exciton operator, 


(j 


t _ 


m, 




(S95) 


np ms 


such that its corresponding transition dipole is superradiandy enhanced at p„ - \/NnpP„. 

Having addressed the effective interaction between a single nanopillar and the SP modes, we can move on to the de¬ 
scription of the superlattice, Hexc-SP - '^nH^J^'c-SP' Pn - P assume periodic boundary conditions (PBCs), 

we can construct Fourier k modes for the excitons too, cr! = ^ Then, we arrive at the Hamiltonian, 

^ \/NxNy 

Hexc-SP = LkHexc-SP.k = Zfc ,/(fc) AfeO-^-t h.c., where 



g-ao(fc)z(fc) 


p-E[k), 


(S96) 


and k runs for all the allowed discretized wavevectors ki---^ + qt for q/ = 0,1, • ■ ■, Ai/ - 1. 

Within the MFA, we have achieved to represent each nanopillar as a single collective transition dipole /r„ associated with 
the operator crj,. There is, however, an ambiguity in this approximation, namely, the criterion to optimize the parameter 
z[k). We will discuss this in Subsec. HUB 31 . 
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2. Beyond the MFA 

Going back to Eq. jS94al and assuming we can consider alternative nanopillar modes. We foUow Gonzalez- 

Tudela, etal . Let us define the new modes (distinguished from the others by the overbar notation), 

aiik) = (S97) 

V>yV ms 

where the corresponding normalization is given by, 


ms 

rzf 


^np 


dz\p^-E{k,z)\ e 


2 -2Stao(fc)2 


(S98) 


where we identified - zo = 144 as the vertical thickness of each nanopillar. Introducing the corresponding k mode = 
Eq. )S94a) becomes Hexc-SP = ~ Zt </(*:)+ h.c., where, 


where we have identified 


Jlk): 


2Stao(fe)z|p^.£(^^_2)|2^ 


(S99) 


NxNyNnp {NxWx){NyWy) 

SM4 “ S ^ ’ 

^ ^- • 

<1 

as the average density of chromophores in the organic superlattice which, due to the “void space” between nanopillars, is 
lower than p„p. Except for a different convention in the phases of our exciton modes this solution has the same structure 
as the one presented in , even though the latter deals with an organic layer of uniform density. 


3. Comparison 

When we wrote Hgxc-SP ~ Zt a?(fc)afccr^ + h.c., we made a MFA to Eq. l lS94al by invoking definitions for the mode crj, and 
the coupling ^{k] (Eqs. jS95l and IS96H . The essence of this approximation is the exponential factor e-“o(fc)z(fc) ^[k], 
which implies that when each nanopillar interacts with the fc-th electromagnetic mode, it behaves as a collective dipole 
placed at the effective height z = z(fc). On the other hand, when going beyond the MFA, we introduced d], (fc) and ^{k] (Eqs. 
IS97> and I IS99H and showed that Hexc-SP ~ Zfca?(fc) + h.c. is an exact representation of Eq. IS94a> (notwithstanding 
the excellent approximations of converting the sums over chromophores to integrals and the PBCs). Hence, d^ (and not cr^) 
is the natural exciton mode which couples to the SP mode 

The solution beyond MFA might be a more convenient description if one is interested in a careful description of the en¬ 
ergetics of the problem. However, for purposes of the topological characterization of the plexciton bandstructure, it is more 
pertinent to adopt the MFA description. Notice from Eq. <5971 that djj(fc) depends explicitly on k, so that the Fourier modes 
d! = Z_ Znh'n(fc)g'*^''^ have an additional dependence on k beyond the phase factor This introduces a technical- 

ity for the numerical computation of the Chern number for each plexciton band, which we wish to avoid at present. This 

4 - J- 

complication does not occur in the MFA, where (t„ uniformly sums over the excitons operators bl^s for a given nanopillar 
regardless of k. Thus, we make a compromise; we formally use the structure of the MFA, but heuristicaUy make the energetic 
approximation \^{k) - ^{k). Gomparing Eqs. IS96> and <5991 . this identity requires that z(k) satisfies. 


p„-£'(fc,z(fc))e"““'-«‘>®^® 



dz|p„-£(fc,z)pe 2Staorgo(fc)z_ 


(SlOl) 
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This constraint has the appeaiing physicai content of computing the mean-fieid effective position z[k) of the coUective 
nanopiiiar dipoie by averaging interaction of the SP with respect to the intervai z e [Zi,Zf]. Operationaiiy, however, it is not 
necessariiy to soive for z{k), as having the vaiue of ^(fc) in Eq. | |S99> suffices for our caicuiations. 

To remind ourseives, when deaiing with the dispersion of the organic iayer alone, we have coarse grained it into an effec¬ 
tive rectanguiar superiattice. However, when computing the exciton-SP coupiing via Eq. IS99i , p must be taken to be the 
density of the originai monociinic superiattice. Eq. l lSlOOi is undetermined as Nx and Ny are artificiai simuiation parameters, 
instead, we shaii write 


p^fpnp, (S102) 

the average density is equai to the surface coverage fraction times the originai nanopiiiar density. 


4. Representative coupling values 

Figs. [S6] and [ST] piot representative exciton-SP coupiing vaiues ^{k) using Eq. iS99l for the two-iayer and three-iayer 
setups, respectiveiy. We dispiay caicuiations for different orientations of the transition dipoies of the nanopiiiars, when 
^l\\k,^l II 0k, and fi II z. Throughout the piots, we have chosen siiver Drude parameters Coo ~ 4, wp ~ 9 eV and g = 0.3. Each 
of the paneis dispiays resuits corresponding to a particuiar dieiectric permittivity and the base height of the nanopiiiars 
Zo, and fixing the same organic iayer height 144 = 70 nm. Taking the parameters for the organic superiattice described in 
Subsec. and assuming pnp - 37chromophores/nm^, the average density using Eq. jS102i is p = 14chromophores/nm^. Fig. 
[SOjcorresponds to the physicai scenario where the organic nanopiiiars are “embedded” in the MO dieiectric iayer starting at 
the base height zq, whiie Fig. [STjassumes that the MO dieiectric iayer has a width a and the base height of the nanopiiiars is 
aiso at a. 

As a reminder, in the absence of the MO effect, the eiectric fieid of the fc-th SP mode has no tangential component (see 
Eqs. IS40ai . IS72aH . Thus, the blue curves in the plots (/r || 6) must vanish identically for g - 0. For g ^ 0, these couplings 
scale linearly with g (see Eos. IS41^IS76i . so it is easy to predict these perturbative exciton-SP couplings for other values of 
g. On the other hand, the red (/r || k) and black (/u || z) curves are independent of the value of g, as they are 0(g°). Notice 
that all the curves peak at some short wavevector kmax- The effective site energy (ii)eff (see Eq. <S87H was optimized so that 
the exciton dispersion Hgxc,k and the SP dispersion Hspk become degenerate (Dirac points) at wavevectors k* satisfying 
|fc* I = kmax- This is carried out so that upon inclusion of the MO effect, the topological anticrossing generated at the Dirac 
point k* becomes as large as possible. The latter implies that one should maximize ^{k*) for fi || Ok- 

The simulations displayed in the main text correspond to the first (upper left corner) panel of Fig. jSS] where, yielding 
^{k*) - 0.24eV for p || 6, or equivalently, a topological anticrossing gap of 2 ^{k*) - 0.48eV. This value becomes reasonably 
large compared to the linewidths of the exciton and SP modes even if we reduce g by a factor of three or four. Note that quite 
often, the largest couplings occur when fi || z, reaching values which are comparable to the exciton site energies. This regime, 
known as ultra-strong coupling |^, is interesting in its own right and gives rise to novel effects, which are beyond the scope 
of our work. Unfortunately, for our purposes, we cannot exploit these large couplings, as they do not vanish for any k and 
hence, does not yield Dirac points (see main text). 

A few interesting trends can be obtained from scanning through and zq; some of these results are displayed in Figs. |S6| 
andjSTj First, couplings decrease as zq increases. This is not surprising, as owing to the evanescent nature of the SP fields, the 
coupling should be strongest for the chromophores that are closest to the interface at z = 0. Second, couplings decrease as 
the dielectric increases. Finally, notice that the three-layer setup yields very weak values of ^{k*) for p || Ok- We believe 
that these two problems are related to index mismatch between the various interfaces. As mentioned in the main text, one 
possibility to ameliorate this problem is to embed the MO material inside of a low dielectric polymer. 

An alternative solution to increase ^{k*) for fi || 0k is to use materials with large g values at the UV/visible, which is what 
we need (recall the crossing between the SP and exciton dispersions at k* happens at 3.1 eV in our calculation). Some exam¬ 
ples of the latter are Co alloy films El , orthoferrites , or spinels . A caveat about the latter is that they are also highly 
absorptive at those same wavelengths ( larg e imaginary part of ea, which we have neglected in this work). Ce substituted 
YIG has less of a problem in that regard We are currently addressing all these possibilities, including different stacking 
geometries, in order to induce strong MO effects. 
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Figure S6. Representative exciton-SP coupling values for two-layer (metal-MO dielectric) setup. The calculations have been carried out 
using Eq. 15991 . We display results for different orientations of the transition dipoles of the nanopillars, when fi \\ k, fi || 0, and /i || z. 
Notice that all the curves have maxima at short wavevectors. The calculations assume that e^o ~ 4, wp ~ 9 eV g = 0.3, and the height of the 
nanopillars being = 70 nm. 
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Figure 57. Representative exciton-SP coupling values for three-layer (metal-MO dielectric-organic) setup. Just as with Fig. the calcula¬ 
tions have been carried out using Eq. (599]. We display results for different orientations of the transition dipoles of the nanopillars, when 
p\\ k, p. II 0, and p || z. Notice that all the curves have maxima at short wavevectors. The calculations assume that Coo ~ 4, wp ~ 9 eV, 
g = 0.3, and the height of the nanopillars being Wz = 70 nm. 
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